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SUMMARY 



The centered difference leapfrog technique, commonly used in the numerical 



solution of the primitive equations of meteorology, computes the gradient 



of the geopotential at the n~^ time step by: 

f£ (t n'V ' 2S { *j+l ‘ 



♦J-l ' • 



where: 



♦j " * (t n- x j ) 



The Shuman technique replaces this term by a weighted average of the centered 

s t th s t 

differences on the (n-1) , n , and (n+1) time steps. This should 

allow the use of a longer time step than that predicted by the CFL condition 
before the onset of computational instability. This paper considers the 
computational stability regions obtainable by this technique in the cases 
of: (1) the linearized barotropic model with no Coriolis force and no mean 

flow, (2) the linearized barotropic model with no Coriolis force, but includ- 
ing a mean flow, (3) the linearized barotropic model with no Coriolis force 
and no mean flow, but incorporating the time filtering technique designed 
by Robert, and, (4) the linearized barotropic model with no Coriolis force, 
but including a mean flow and time filtering. In each of these cases, the 
Shuman technique is demonstrated to yield an increased time step and the 
maximum size of the time steps is shown. The presence of a mean flow or 
filtering yields maximum time steps lower than those of basic solution 
(case (1)), but still significantly higher than the unaveraged method. 
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1. Introduction 



Shuman [6] has proposed a modification to the leapfrog scheme for the 
primitive equations which are used in numerical weather prediction* The 
scheme which uses a weighted average of the pressure gradient at the past, 
present and future times, allows a longer time step than the one given by 
the CFL condition. Brown and Campana [2] have carried out the linear 
analysis of this modification for the three cases: (1) the barotropic model, 
(2) a two-layer model whose pressure is the vertical coordinate, and (3) a 
two-layer model with the Phillips [4] sigma coordinate system. In this 
paper we will treat only the linearized barotropic model. We will present 
an analytic solution for the computational stability curve for the simplest 
case. A mean flow will be added to the equations and the stability curve 
will be obtained numerically and also with an approximate analytic pro- 
cedure. Finally the time filter designed by Robert [5] will be applied to 
all variables and the resulting stability curve will be presented. It 
will be seen that both the presence of a mean flow and the time filtering 
reduce the time step, but that the use of the Shuman technique is still 
advantageous . 

2. Basic solution 

The linearized equations for the barotropic model with no Coriolis 
force and no mean flow are: 



Here u is the velocity component in the x-direction and (p/g is the 
departure of the free surface height from its mean value $/g . These 




( 1 ) 




( 2 ) 
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equations describe shallow water waves which move with speed . 

The finite difference approximations to equations (1) and (2) with the 
Shuman pressure gradient averaging included are: 



n+1 n- 
u . -u . 

_J L 

2 At 






,n+l m- 1 



C- g 

2At 



n n 

u.,, - u. . 

( _J±L lzl) = o 

v 2 Ax ' 



(3) 

(4) 



where the discretization uses x = jAx and t = nAt . The usual leapfrog 
differencing is obtained by setting a = 0 . This scheme is explicit since 
may be obtained from (4) before it is needed in (3). 

In order to obtain the computational stability properties of this 
scheme we substitute the following expressions 



n . n 
u . = A co 
J 



ikAxj 

e 



(5) 



j.n n ikAxj 

<p. = B co e 

T J 

into (3) and (4). After the constants A and B have been eliminated we 
obtain the following quart ic equation for co: 

co 4 + 4 S a co 3 -2[l-2S(l-2 a )] co 2 + 4S a co+1 = 0 , ( 6 ) 

where S h ( 7^) 2 $ sin 2 kAx . (7) 



Since S and a are real, the roots of ( 6 ) are either real or in complex 
conjugate pairs. Due to the symmetry in the differencing we expect that 
|co | = 1 for some range of the parameters. As a result the polynomial 
factors to 

(co + 2a^a) + 1) (a> + 1) = 0 . (8) 
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If we expand (8) and compare with (6) we find that a^ and a must 
satisfy the following equation: 

a ± 2 - 2Sa a ± + S(l-2 a ) - 1 = 0 , i = 1,2 . (9) 

The solutions to this equation are 

= Sq +|/ (Sa+1) 2 " S . (10) 

It can be seen from (8) that only when the a^'s are real and of magnitude 
less than or equal to unity will all the roots of (8) have jco | = 1 . The 
value of a^ will be real when the quantity under the radical is non-nega- 
tive; therefore, this condition can be written 

(S a + l) 2 - S > 0 . (11) 



If we choose the equality we will obtain the following stability relation: 



S = 



1 



2q - a A - 4q 







( 12 ) 



where the minus sign in front of the radical was chosen to make the expres- 
sion reasonable in the limit as q -* 0 . The condition that the a^'s have 
magnitude not greater than 1 leads to the condition 




(13) 



This is consistent with (12) which becomes complex for q ^ • 

The combination of conditions (12) and (13) is given in Fig. 1 as the 
curve labeled o - 0 . Brown and Campana [2] determined a curve with 
approximately this shape by numerically solving for the roots of (6); 
however, they did not present its analytic form (12). The maximum value 
of S which is S = 4 , occurs at q = ^ . For the usual leapfrog 
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differencing (a = 0) the maximum value of S which allows computationally 

2 

stable solutions is S = 1 • Since S is proportional to At (see Eq.(7)) 
it follows that the use of the Shuman pressure gradient averaging in this 
linear system allows a doubling of the time step as compared with the 
standard leapfrog scheme. However it may be difficult to achieve this 
factor of 2 in practice when other effects are included since the width of 
the stable region goes to zero as S approaches 4 . In fact a value of 
a which is slightly less than ^ would probably be preferable. 

3. Solutions with mean flow 

In this section we add the effects of a constant mean flow to the 
stability analysis. The linearized equations are: 



Su 

ST 



+ u 



Su . 
5x 



dx * 



(14) 





(15) 



Normally the addition of the mean flow terms does not have a great effect 

on the computational stability criteria since the phase speed of the 

1/2 

external gravity waves, $ , is generally much greater than the speed 

of the mean wind, U . 

When Eqs. (14) and (15) are put in finite difference form with the use 
of the approximations in (3) and (4), the mean advection terms are evaluated 
with centered finite differences, and the relations (5) are substituted 
into the finite difference forms of (14) and (15), we obtain the following 
equation for to : 

cd 4 + 4(Sa+ la) 2[2S(l-2 a )-(l+2a 2 )] u> 2 + 4(S a - ia) a> + 1 = 0 , (16) 



where 



UAt 

a = 



s in kAx . 



(17) 
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Note that dividing (16) by co and taking the complex conjugate of the 
equation, yields an equation in 1/co* with the same coefficients as (16). 
Thus, it follows that a) is a root of (16) if and only if (co*) ^ is 
also a root, and therefore the roots must be distributed in exactly one of 
the following ways: 



Case I: 



Case II: 



Case III: 



Four roots (including multiple roots) on the unit circle 
(e^l, e 1 ^, e 1 ^, e X ^4; + t 2 + $3 = ” * 

Two roots on the unit circle, two roots off the unit circle 
(e 1 ^, e^3, pe^l, £ e^l; 2<^ + ^ = - <^) . 

Four roots off the unit circle 
(pe 1 ^, ^ e 1 ^, £e""^, j e" 1 ^, p,£ i 1) 



In each case the fact that the product of the roots is equal to the last 
term in (16) was used. 

The analysis of Cases I - III indicates that (16) has at least one 
factorization of the form 



, 2 , _ i0/2 

(co + 2ae co 



+ e 10 ) 



t 2 . -ie/2 

(o) + 2be a) 



+ e" i0 ) = 0 



(18) 



for some angle 0 , and a and b real. Furthermore it is easily shown, by 
analyzing the possible combinations, that in either of the unstable cases 
(Cases II or III), the factorization (18) must be unique and either a or b 
(or both) greater than unity in magnitude, whereas the stable case will 
generally have three distinct representations with both a and b less than 
or equal to unity in magnitude. The only exception will be when stable, 
but degenerate (multiple) roots occur. 
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If we expand (18) and compare with (16) we can obtain the following 
equations : 

a - Set sec 0/2 4- q esc 0/2 , (19) 

b = Sa sec 0/2 - a esc 0/2 , (20) 

S 2 a 2 sec 2 0/2 - S + s a + cos 2 0/2 - o cot 2 0/2=0. (21 ) 

2 

When we write u = cos 0/2 and use trigonometric identities we can obtain 
the following polynomial in u : 

P (u) s (u-1) (u 2 + S(2a-1) U + S 2 a 2 ) + a 2 u 2 = 0 , (22) 

a 

where one different factorization of the form (18) arises from each distinct 
solution of (22) in 0 < u < 1 . Thus by our above comments a sufficient 
condition for stability is that (22) have more than one solution in [0,1]. 
The appearance of (22) as a cubic should be expected, since as noted above 
in the generalized case of stability, we expect three distinct factoriza- 
tions. When only one solution of (22) occurs in [0,1], then there must 
be instability unless (16) has a triple root. It can be shown that there 
is at most one value of q for which (16) has a triple root. 

With this characterization consider the qualitative behavior of P (u) . 

a 

Note for all u > 0 , a > 0 that P (u) will lie above the curve 



' a 

P o (u) = (u-1) (U 2 + S(2a-1) U + S 2 a 2 ) , (23) 

where it is easily seen that 

p o (0) = P o (0) = - S 2 a 2 , (24) 

P o (l) = 0 , (25) 

P' (1) = (Sa + l) 2 - S . (26) 

o 
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Also note that for fixed u, P (u) is a strictly increasing function of a • 

o 

These observations now allow us to visualize fairly easily how an in- 
stability develops. Consider Fig. 2, which shows a typical situation when 
P q (u) has three distinct roots in 0 < u < 1 , i.e. o - 0 yields a 
stable procedure. Then observe that as a increases the curve moves 
upward, causing the value of the smaller root to decrease, and the two 
larger roots to move closer together. Eventually, for some value o , 



the larger two roots degenerate to a single double root with P 



(u) 



max 

tangent to the axis there. Finally for a > c mav , the curve detaches 
from the axis, leaving only the single smaller root, and hence instability 
must occur. 

Intriguingly, this construction also allows us to visualize a situa- 
tion when a > 0 will yield a stable procedure even though a = 0 does 
not. This is shown in Fig. 3, where as a increases the curve first 
attaches itself producing a double lower root, which then moves apart, 
until eventually a double upper root appears, followed by instability. 
Numerical solution of (16) verifies this in fact occurs, as evidenced by 



the curves for a > £ in Fig. 1, 



Now consider the situation when P^(l) > 0 and 0 < a < ^ 



of P q (u) = 0 are 



U o = 1 ’ 



r<l-2a) - V 1 - 4 a ] , 



u 2 - f [(l-2a) + Jl - 4a ] 



The roots 

(27) 

(28) 
(29) 



Since (l-4a) < (l-2a) it follows that P q (u) has a shape as indicated 

in Fig. 2. The maximum stable value of a (hereafter called a ) occurs 

max 

when the larger root of P' (u) =0 is also a double root of P (u) . 

a o 
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Since P^(u) a cu bic the conditions can be worked our formally; how- 

ever, the resulting conditions are too cumbersome to provide computational 
or analytic insights. Therefore we shall present some simpler necessary 
conditions and approximate expressions. 

First, observe that if P^(u) ^ 0 > P Q (u) has exact ly one root in 

0 < u <1 , hence P (u) must have exactly one root in 0 < u < 1 . Thus 

a 

P^(u) >0 is a necessary condition for stability when a > 0 • (Observe 

the condition that (26) be non-negative is identical to (11).) 

A rather accurate approximation to a for this situation can be 

max 

obtained by noting that the onset of instability must correspond to a = 1 , 

since a > 1 yields, from (18), a solution with |cu| > 1 • Thus from (19) 

a satisfies 
max 



or 



Sa sec 0/2 + a esc 0/2 = 1 , 
max y 

= 1 , (30) 

A ^" U L 



where u T denotes the smallest root of P (u) . This is necessary 
i* o 

max 

since both a and b must be continuously dependent on a . 

The u in (30) is not known exactly, but the examination of Fig. 2 

Li 

indicates that 

U L ~ u i > (31) 



where is the lower root of P q (u) ♦ Thus to a first approximation, 

stability will occur for 



0 < a = J 1 - 
max v 



u. 




0 



1 

< a < ^ > 



(32) 
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where is given by (28). It will be seen that this equation gives an 

excellent approximation to the values obtained numerically. Unfortunately, 
a similar approximation for ^ < a has not yet been found. 

The roots of (16) were computed numerically and the resulting stability 
curves are given in Fig. 1 for selected values of o • All the curves in 
Fig. 1 have similar shapes with the lowest curves corresponding to the 
highest values of a • The following stability condition for a = 0 is 
easily obtained: 

S 1/2 + |o| < 1 . (33) 

For a = 0.1 the increase in At over the value for a = 0 (see equation 
(17)) is about 80 percent. A typical value of q for operational numerical 
prediction models would be between 0.1 and 0.2. Fig. 1 also shows that 
stability can occur for the larger values of a in the region a > £ • 
However these values of a should not be used because there will always 
be a value of kAx in (17) which will give an arbitrarily small value of 
G and therefore instability. 

Tables 1 and 2 compare values on the numerical stability curve from 
Fig. 1 with values computed from (30). The agreement is quite good with 
the largest difference occurring for a near ^ and for larger values 
of Q 

4. Solutions with time filter 

Time filtering has been used in numerical weather prediction models 
to damp both physical and numerical noise (Robert [5], Haltiner and 
McCollough [3]). Consider the following centered time filter: 

F(t) = F(t ) + 7 [ F(t + At) + F(t - At) - 2F(t)] . (34) 
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This form is convenient for operational prediction because it uses the 
previous averaged value which saves machine storage. When this filter 
is used in linear equations which have solutions proportional to a / 1 , 
equation (34) takes the form 



F(t) = (F (t ) +r[F(t + At) - 2F(t)]/(l - roT 1 ) 



(35). 



Here we have used the relation F(t - At) = 03 1 F(t) . 

This time filter is introduced into the finite difference equations 
(3) and (4) by replacing u™ ^ and (j)^ ^ with the filtered values ob- 
tained from (35). When the relations (5) are introduced into the time 
averaged difference equations we obtain the following equation: 

2 



H 



to 4 + 4 [S a -r] aj 3 + 1 4 ( 7(1+7) + s[l-2a(l+7)] 



)- 2 ] 



03 



+ 4 [S (a ( l+2y) “2/ ) + 7 ( 1 - 27 )] CD + 4S7(7-a) + (1-27) = 0 . (36) 



If we consider the special case of no pressure gradient averaging 
(a - 0), the equation reduces to 

[co 2 - 27 cd - (1-27 ) ] 2 = -4S(cd-7) 2 . (37) 

Take the square root of both sides of this equation and solve for 03 
which yields 

CD = 7 ± iS 1/2 + y (7-1) 2 - S . (38) 

This result was obtained by Asselin [1] who has discussed the solutions 
in detail. When S is sufficiently small the solutions will be damped. 
The critical value of S which is always less than 1, decreases as 7 
increases. However, in the damping region, the damping rate increases 
with increasing 7 . 
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In the general case the roots to (36) must be found numerically. 

Fig. 4 contains the curves which separate the unstable solutions from the 
stable solutions for selected values of 7 . The left hand limits of the 
curves show the reduction in the critical S as a function of 7 for 
a = 0 . The curve for 7 = .05 closely approximates the curve for 
q = 0 in Fig. 1. As 7 increases the maximum stable value decreases and 
shifts to the right. In fact sizable stable regions exist for a > ^ 

depending on 7 . For 7=0 there are no stable solutions for a > ^ • 

We now consider the effect of the time averaging on the solutions when 
the mean flow is included. When the time averaging effects are added to 
equation (16) we obtain: 

cjdS- 4 [Sa-7 + ia]u^+[4(7(7+l)-a^+ S (l-2a(7+l) ) )-2 -12a7i] to 2 
+ 4 [S(a(l+27)-27)47 (1-27+27 2 ) + ia(27(7+l)-l] co 

+ 4S7(7(l-CT 2 )-a) + (1-27 ) 2 + 14(37(1-27) = 0 . (39) 

Fig. 5 contains the stability curves which were obtained by numerical solu- 

tion of (39) for the value a = 0.1 . The curves which are for 
7 = 0.1, 0.2, 0.3 show smaller values of maximum values of S than are 
seen in Fig. 3. However the starting points (a = 0) are also smaller. 

In general the peaks are located at about the same values of a and the 
peaks are broader. Both Figs. 4 and 5 show a large improvement in maximum 
S over the value with no pressure gradient averaging (a = 0). 

5. Conclusions 

The stability properties of the Shuman [ 6 ] pressure gradient averaging 
technique have been investigated in this paper with the linearized shallow 
water equations. The analytic solution obtained in Section 2 shows that 
the time step can be doubled for a = jr • However, the width of the stable 
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region becomes very narrow as a = ^ is approached so that the best value 
of a would be slightly less than ^ . When a mean flow is included the 
time step must be reduced. However for reasonable values of the mean flow 
(q = 0.1 to 0.2) the time step can still be increased by 70 to 80 percent. 

The time averaging of all variables which was suggested by Robert [5] has 
been used to damp unwanted high frequency components in numerical forecasts. 
The use of the time filtering, however, requires a smaller time step and 
therefore more computation time. When the time filtering is used in con- 
junction with the pressure gradient averaging, the time step can be 
significantly increased although for the larger values of 7 the time 
step may not be much larger than with no time or pressure gradient averag- 
ing. When the time averaging is used the optimum value of a is critically 
dependent on 7 . The addition of the mean flow decreases the time step, 
but does not appreciably affect the optimal a • 

The Shuman [ 6 ] pressure gradient averaging technique has been used 
operationally at the National Meteorological Center and it is now under- 
going tests at the Fleet Numerical Weather Central. This technique should 
be useful in other fluid dynamical applications provided that the 
velocities are appreciably less than the fastest gravity waves. 
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TABLES 



a 


0 


0.1 


0.2 


0.225 


0.25 


s 


0.78 


1.00 


1.50 


1.75 


2.62 


an 










S 


0.80 


1.02 


1.55 


1.80 


2.50 



num 



Table 1. Comparison of the values of S on the stability curve from 

Fig. 1 (S ) with the values of S computed from Eq . (30) 

(S ) for n a m = 0.1. 
an 



a 


0 


0.1 


0.2 


0.225 


0.25 


s 

an 


0.33 


0.43 


0.61 


0.71 


1.04 


S 

num 


0.35 


0.46 


0.62 


0.74 


0.88 


Table 2. Same 


as Table 1 except 


for a - 0.4 . 
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20 
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